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Abstract 


In the; present work, a tliree-clinrensional large deformation elasto-plastic finite element 
code has heem ckweloped for the analysis of deep drawing process. The parametric study 
is in'cseiited t.o show the effects of two process parameters namely sheet thickness and 
ma.teria.l propc'rties. 

The up(UU.('(l Lagrangian formulation is used to develop the finite element equations 
using H-n()d('d solid eh'iiK'iit. Tlu' new incremental objective stress measure' and loga.- 
ril.hmic sl.ra.in measure are used which allow the use of large incremental displacement. 
Modific'd Ne'wton-Raphson iterative technique is used to solve the nonlinear incremen- 
tal e(iua.tions. Tlu' material is assumed to be elastic-plastic strain hardening yielding 
a.ccording to l.lu' von Mises criterion. The strain hardening behavior is modeled by a 
I)OW('r la.w. Diu' to small a.ccelerations, inertial forces are not included. The body forces 
arc,' also nc'giected. The friction is assumed to be of sticking type at the punch-sheet 
interlVici'. while the die-sheet and the blankholder-sheet interfaces are assumed to have 
.sliding friction. To reduce the storage and computational time, the skyline method of 
assembly a,nd the static condensation of coefficient matrix are carried out. 

Va.iida.tion of tlu.; code is done by comparing the predicted finite element results with 
a,va.ila,bl(' expi'rimental results. Parametric study shows that the punch force increases 
wil.h t.he increa.se in slu'ct thickness. The thickness strain is obeserved to be the highest 
in t.he clea.ranee region while the equivalent stress is found to be the largest at the punch 


corner zoiu'. 
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Chapter 1 


Introduction 


Sheet metal parts are one of the predominant commodity in different engineering ap- 
plications. Manufactnring of sheet metal parts by means of press working is a cost 
effcHdive procicss since it eliminates expensive machining and welding operations giving 
better quality finished product and it enables to produce components at a very high rate. 
Among the manufacturing of these sheet metal parts, deep drawing is an extensively used 
press working process. 
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Figure 1.1: Deep Drawing Process 


De(‘j) drawing is a process in which a blank or work piece, usually controlled by 
a pressure i)latG, is forced into and/or through a die by means of a punch to form a 
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hollow compoiieiii [1]. The procedure with greatest range of application is the deep 
drawing witT rigid tooling, involving draw punches, a blankholder and a female die. 
Tlie blank is i)ulle(l over the draw punch into the die, the blankholder prevents any 
wrinkling taking place in the flange [2]. The schematic diagram of the process is shown 
in figure (1.1). Deep drawing without a blankholder is carried out only if the part does 
not tend to wrinkle, i.e. it is thick or has small amount of deformation [3]. The workpiece 
is subjected to radial tension forces and tangential compressive forces. Therefore, the 
material is (‘oiniu’essed in the tangential direction and stretched in the radial direction. 

Dcei) drawdng of non circular components finds applications in many industries, 
such a.s construction, automobile manufacturing, aerospace technology, and electrical 
equipments. However, the stress and strain behaviors exhibited in non-circular deep 
drawing j)ro('(‘sses are so complicated that the mechanics of deformation differs very 
a.i:)i)re(‘ial)ly from that ol)served in circular deep drawing process. Among the non circular 
dcK']) drawing pi'occ^sses, box drawing is the most basic process. A major problem in box 
drawing is t.ht' diffca'cnce of metal flow rates between the straight walls and the corners 
of t.lie deforming j)art. This will result in an uneven material distribution around the 
Clip wall. In axldition, the variable flow rates may also cause wrinkles in the flange of 
the drawn ciup. 


1.1 Literature Review 

The development of reliable analytical procedures to predict the deformation and stress 
fields in shcct-metal forming processes is a possible solution to make the process more 
efficient. However, this development has encountered serious obstacles. The major cause 
is t.he i)r(,;sc;iice of various non-linearities such as 

• the geoim!l.ric non-linearity due to large deformation of the readily pliable sheet, 

• the' ma.t(;rial non-linearity due to the deformation going into plastic range, and 

• the kinematic non-linearity due to the contact of sheet metal with punch and die 
along wit,h friction at the contacting surfaces. 

K4(>d('liug nf tlu'se non-linearities is a difficult task due to the complexities involved in the 
midciiying mathematics and the numerical instabilities that come up during the solution 
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stage. 


Ill the jia, str, a iiumbei' of approximate methods of analysis viz., the slab method, the 
slip line field method, the visio-plasticity method, the upper and lower bound techniques 
and the Hill s geneial method, have been developed for various forming processes [4]. 
In the slab method, the workpiece being deformed is divided into several slabs. For 
each slab, simplifying assumptions are made with respect to stress distributions. The 
resulting approximate equilibrium equations are solved with imposition of stress com- 
patibility between slabs and boundary tractions. The solution of these equations gives 
approximate load requirement and stress distribution in the forming process. The slip- 
line field method is used in plane strain problems for perfect plastic materials and uses 
the hyperbolic properties of the stress equations. The construction of slip-line fields, 
although i)ro(hK:eH an “exact” stress distribution, has limitation in predicting results 
compa.rable with experimental work. From the stress distributions, velocity fields can 
be c:a.lcula.t('d tJu'oug'h i)lasticity equations. The visio-plasticity method combines experi- 
mc'iital a, ml a,na.lytical methods. After the velocity vectors are determined from an actual 

sti'a.iu~ra,t(',s a, re calculated and the stress distributions are obtained form plasticity 
(K|uations. Th(' method has helped to obtain reliable solutions for the processes where the 
cx;)erimenta.l dc'tfu'inination of velocity vectors was possible. The upper-bound method 
retpiires the giiessing of admissible velocity fields, among which the best one is chosen 
by minimi'/ing total energy. Information leading to a good selection of velocity fields 
come from exix'rimcntal evidence and experience. This method, with experience, can 
dediver a fast and relatively accurate prediction of loads and velocity distributions. Hill 
has given a gtmeral method of analysis for metal-working processes when the plastic flow 
is unc;onstrain(ul. The method is based on a criterion of approximation derived form the 
intcn-pnjtation of the virtual work-rate principle. 

I'lu' a, hove irndbocls have been useful in predicting the forming loads and overall 
gcoiiK'try chauj!,cK of deforming workpieces with approximation. However, accurate de- 
termimUaon of t!>c effects of various process parameters on the detailed metal flow has 
Ixicome possible with the use of finite element methods. This method is sufficiently 
versa, t,ile to ha.ndl(i various geometries, materials and boundary conditions, hence widely 
used in t.he a,nalysis of sheet metal forming. 
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1.1.1 Finite Element Method 


In the simviLU.ion of sheet metal forming, several finite element models such as membrane, 
shell and solid (3-D) models have been proposed. The membrane model is the simplest 
one. Toll and Kobayashi [5] modeled the square cup drawing process based on finite 
stiain foiniulation and membrane theory. An optimum blank shape was determined 
based on finite element lesults of net material flow. But they assumed zero blankholder 
force in the computations. Saran and Samulesson [6] modeled the behavior of sheet 
materials by a t.riangular membrane constant strain element. They used hypoelastic- 
viH(!oi)la8(.ic matia'ial model including Hill’s plastic anisotropy, power law hardening and 
strain ra,te smisitivity. They demonstrated that, a huge reduction in the number of 
eleuK'nl.s can bc' achieved by adaptive re-meshing. Majlessi and Lee [7] developed a 
finite (’lenient a.nalysiH technique using elements with membrane characteristics. They 
ignored the' (‘ffects of bending and unbending and found approximate solution to the 
problem of noiuixisymmctric deep drawing. However, the membrane formulation neglects 
,stre.s.s variation in thickness and considers the stress components only in the plane of 
sheet, .'I'hus, il, i.s not ajipropriate to erhploy the membrane theory for conditions in which 
bending deformation of sheets becomes significant and dominant like in deep drawing. 
In general, the iirocesses where the bending effect is significant are simulated by the shell 
mocU'l or tlu' solid model. 

Onatc and Saraciber [8] presented the viscous voided shell approach for developing 
the bending/membrane finite element formulation for sheet metal forming. These viscous 
voided shell elements can be selectively used for membrane or full bending analysis of 
the slu^cit according to the nature of deformation. They concluded that for 3D cases, 
using compk^x ptmehes and dies, this approach is effective and economical. Chou et 
al. [9j sinmlat.ed the sheet metal forming by plane-strain shell elements using the stress 
resultant l,heory. This theory is based on quadratic yield function, a hardening rule, 
and (,ht> a.ssocia,1.(xl flow rule. Shi et al. [10] introduced a new finite element approach 
for the diixict i)rediction of blank shapes and strain distributions for desired final shapes 
in slKct niet.al forming, based on a one-step simulation algorithm. They have given 
new a,nd c:om])act formulations of the DKT^ shell element. Although the shell model is 
more (dibctivc’ than the solid model, it still takes a substantial amount of computational 
iinui and imnuory spacie for its integration in the thickness direction. Chou et al. [12] 
propos(‘d il, sl,reHs resultant constitutive law based on the orthotropic yield function of 

* '''‘ lVis<’ret('’ 'Kir(li<')(rTla,u^ (DKT) shell element, is widely employed in non-linear static and dy- 
iiiunic analysis [11]. 
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Hill fot slieotiK with plaiiai isotropy. They avoided integration through-the-thickness by 
coiisi(l(!i iug ( 1 . i)la.to element, and showed a substantial amount of saving in computational 
time. 


Keck et, al. [13] modeled the deep drawing of rectangular blanks using eight noded 
hexahedials. They assumed a constant gap blankholder force in the analysis, but in 
practice the l,)laiikholder force is given by the die cushion. Lei et al. [14] modeled an 
elasto-plastic finite element solution based on solid element and finite-strain plasticity 
for a {iv.c bending and a square cup drawing process. The law of constant friction is 
assumed a.t t.lu^ die-workpiece interface and the coefficient of friction remains constant 
during the ])rocc!ss. They investigated the punchload, metal flow, stress distribution and 
Hi)ringl)a,ck (iua.ntit.atively. They have shown that for small initial workpieces, it is better 
to us(' a, circuhu' l)lank shape , and for large initial workpieces, an optimum blank shape 
is mor(’ d(!sira.bl('. Kaiping et al. [15] simulated the square-cup-deep drawing using an 
8-nod(' 3-1) hi'ick edement with one integration point which is a mixed element based on 
the’ Hu-Washi'/,u principle with hourglass control, and a 4-node quadrangular 3-D shallow 
shell (dcuKail, with six degrees of freedom per node [16]. They showed that for the same 
mc'sh a,nd the sa,mc simulation the 8-node 3-D mixed elements have consummated much 
less (,'PU time’ a.nel usexl Iciss loading steps than the 4-node 3-D shallow shell elements. 
Thus in some’ H.i)plications, the shell elements might be replaced by 8-node 3-D mixed 
edennemts without, leasing much precision and with more economically. Explicit finite el- 
(unent teduiiepie' is a dynamic procedure employing the central difference method. The 
advantage of t.his technique is that, there is no need to solve the system of equations and 
therefore, considerable smaller CPU time is needed, as compared to implicit technique. 
Mamalis et a.l. [17] simulated the square cup deep drawing process using the explicit 
uon-liuefi,r finit(’ element code DYNA-3D. Error estimation and mesh enrichment proce- 
dur(!S a.rc' u{’C('ssa.ry to define accurately the current shape of the sheet for the simulation 
of complex forming i)roc(’KS. Bonet [18] showed that these procedures can be efficiently 
imph’inent.cid in existing finite element software with negligible overheads. But, they 
did not, considc’r chirefineraent and removal of nodes from regions exhibiting small errors 
which is (;ss(^ntial for the processes like deep drawing. Fine die details like sharp coi- 
iH’rs (uui not, !)(■ simulatcid with the above techniques. Menezes and Teodosiu [19] have 
publislK’d a, deep drawing simulation of solid model using 8-noded brick elements using 
a,ugmenl,ed La,gra,ngiau method. Their formulation based on Jaumann stress tensor 
a,nd (hven-largi'angt^ strain tensor. They however did not present the variation on punch 
l()a,d a,nd disphuH’nu'ut a,nd also were unable to validate the thickness strain with that 
of t,h(> expr-rimental result. They reported that although larger computational resource 
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was needed tor solid elements the deformation obtained was more realistic. The number 
ol elemcut.H icciuiied to give lesults for solid element were also lesser in comparison to 
shell elernent.s. 


1.1.2 Blankholder Force and Blank Design 


Blaiikholdei force plays an important role in the deep drawing process since it controls 
the material flow into the die cavity. The predominant failure modes in deep drawing 
process arc wrinkling and fracture. In many cases these defects can be eliminated by 
appropriate c.ontrol of blankholder force. Osakada et al. [20] proposed a method for 
determining an oirtimal blankholder force which does not cause wrinkling and localized 
thinning. In this method they used a control algorithm to the FEM program which 
controls t.he l)la.ukholder force depending on wrinkling, thinning, indentation and thick- 
ening. Ahmet.oghi ot al. [21] conducted experiments on an aluminum alloy and arrived at 
bUmkhohU'i' profile a.s a function ot stroke that prevents wrinkling and fracture. Lorenzo 
et al. [22] i)ropo.s('cl the optimal Irlankholder force path which permits to obtain the max- 
imum height, component avoiding both wrinkling and tearing. They have suggested a 
closed-loop eont.rol system based on the fuzzy reasoning interfaced with a FEM code. 
The conf.rol systcmi interacts with the numerical simulation, carries out a continuous 
monitoiing of some relevant process variables like draw-in, punch force, blankholder 
lifting-up and suggests the most effective adjustment of the blankholder force according 
to the imphnneuted knowledge. 


Blank design is a critical part in the deep drawing process. Inverse approach is a 
metliod to deten-mine t.he position of the material points in the initial blank from the 
(!orrc!si)onding m!i.t.eria,l points on the final product. Guo et al. [23] presented the inverse 
a,i)proa.eh to ('Viduat;(! the large jrlastic strains of tridimensional metallic sheets obtained 
by deep dra,wing proc(!HS. Tlu'y used triangular membrane element and deformation 
tlu'oi'y of i)la.st.icit,y. Backward tracing is a key concept to trace backward from the final 
desirable configuva.i.ion to an intermediate preform or initial blank. Ku et al. [24] applied 
t.he ba.ekw'a.rd t,r!i.cing scheme of the rigid-plastic finite element method to the three- 
dinK;nsi()na,l blank dcisign of sheet metal forming. They concluded that sound products 


ca.n be obt.a,ined using backward tracing without any machining after forming. Kim and 
Hull [25] also a.i)plied inverse finite element approach to multistage deep drawing process 
1.0 det.ermine the oi)timum blank shape from the desired final shape. The analysis results 
suf’gesl. t.lia,l. the inl.c^nncdiate die shapes and punch radius should be changed to get the 
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optimum final shape with uniform thickness distribution. Colgan and Monaghan [26] 
tried to determine the most important factors influencing a drawing process, utilizing 
the help of a design of experiments and statistical analysis. The FEA program AutoForm 
was used to model the cup formation. They made the ANOVA^ table which, for a given 
analysis, helps, to determine, which of the factors need control and which do not [27]. 
By their work they have shown that the main aspect that the ANOVA highlighted is 
that the geometry of the tooling is generally most important, especially the die radius. 
But more work need to be done to refine the model and to quantify friction values for 
the punch side and die side of the blank. 


1.1.3 Anisotropy 


Most, shcHit, m('t,a.ls exhibit anisotropic behavior in their deformation. Anisotropy in sheet 
met, a, Is conu's in mostly due to cold working especially during its rolling. The quadratic 
yield crit,eriou by Hill [28] has been the most popular choice to represent the anisotropy 
of sh(Hd,s, particularly for steel. His yield criterion was suitable for rolled sheets, which 
exhibit, cxl orthotropic symmetry. The Hill’s yield criterion and its variants were com- 
fortably applitid to steel and other metals having FCC or BCC type of crystallographic 
str\rcture, but, had anomalous behavior toward highly ductile materials such as alu- 
minum. Yoon et al. [29] cite how this anomalous behavior led many researchers to 
develop other tyi)cis of yield criteria. Hu and Wang [30] describe a selection procedure 
for various anisotropic factors and give an insight on the anisotropic characteristics of 
materials in various metal forming processes. Several non-quadratic yield criteria were 
later developcxl by Hershey [31], Gotoh [32], Banabic et al. [33] and Itskov and Askel [34] 
etc to simulate the anisotropy in sheet metals. In most of these yield criteria the equiva- 
kmt, Htrc\ss w;ls calculated using Cauchy stress and its deviatoric part and anisotropy was 
lai.ei- providfid to the yield function. Each criterion thus developed has its own merits 
for sp<x:ific^ a.pi)lications. 

An (^xt.cmsion of von Mises yield criterion was proposed by Hershey [31] for pol- 
crysta,lLs and lat,er generated by Hosford [35]. This Yield function makes it possible to 
roi)r(!H(;nt yield surface between Mises and Tresca. 

An im|)orta,nt extension of this yield function to anisotropic materials was done by 

'^Th(‘ ANOVA i.s the statistical treatment most commonly applied to the results of the experiment 
to (hiU'nniiu' the pc’rccnt contribution of each factor. 
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Bai'lfi.t et. al. [J(i] which later became popular as Barlat’s yield criterion. They provided 
auLsoti'opy l.o t.hci stress values in different directions and made the yield function an 
intrinsic: luncl.ion of these factors by utilizing the modified stress deviator instead of de- 
viatoric [)ort.ion of Cauchy stress. Karafillis and Boyce [37] later proposed a modification 
to this yield criterion and defined the effective stress to be a weighted sum of two yield 
criteria. Banalric et al. [33] developed a non-quadratic yield criterion for orthotropic 
sheet metals assuming plane-stress conditions. Their yield criterion was a variation of 
Barlat’s. Yoon ct al. [29] presented in their review, a variation of Barlat’s yield criterion 
foi' 2D shec'l. (h'formation. Barlat et al. [38] improved upon this and determined tlie 
anisotropic hectors for a few of the metal alloys in their paper. Bron and Besson [39] 
gave' another la,ct,oi' for anisotropy in Barlat’s yield criterion and also proved for the first 
time the convcxit.y of the yield criterion. They developed a optimization function to 
detf'rminc l.hc iuii,sotro[)ic factors applied to the yield function. They also developed an 
algorithm to imi)lem('nt. the Barlat’s yield criterion to 3D finite element formulation. 

'f’lu' (l('('p drawing ijrocess is one particular sheet forming operation on which several 
r('S('a.rch('i’s hav(' t.ric'd these yiedd functions. The Hill’s yield function has been applied to 
dc’ci) drawing [jrocesH for shell, degenerated solid elements and solid element by Kawka 
a.ud Makiuouchi [•10], Wu [41] Menezes and Teodosiu [19] and others. Barlat’s yield 
criterion has Ix'i'u widely used for deep drawing process especially for shell and membrane 
eU'mcnt.s l:)cca.us(' of its superlative quality in simulating deep drawing in comparison with 
Hill’s yiedd critcilon. Yoon et al. [42], Bron and Besson [39] and Nakamachi [43] are some 
of them. 


1.2 Objective 


The ol)j('ct.ivc of t.hc present work is to develop a 3-D large deformation elasto-plastic 
tinitu eh'ment, code' for the analysis of deep drawing process. The updated Lagrangian 
formula.tiou, tdnd. is convenient for handling material and geometric nonlinearities, is 
used. 'The new incrc'uu'ntal objective stress measure and logarithmic strain measure are 
used instead of idu' .laumann stress rate and Green-Lagrange nonlinear strain measures 
u.s<‘d in most, of the literature. Modified Newton-Raphson iterative technique is used to 
solve (,he uonlinc'ar incremental equations. To minimize the storage, skyline method of 
assembly is used. St.atic condensation scheme is employed to reduce the computational 

tinu'. 
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The material is assumed to be elastic-plastic strain hardening and yielding accord- 
ing to von Mises yield criterion. The effects of temperature and strain rate (viscoplastic- 
ity effects) on the yield strength of the material are ignored in this work. The inclusion 
of these effects renders the analysis quite complex. Due to small accelerations, inertial 
forces are not included. The punch and die are assumed rigid. The friction is assumed 
to be sticking type at the punch-sheet interface, and sliding at the die-sheet interface. 
The blank holding force is applied incrementally assuming to reach the total at the end 
of all increments. 

First, the code is validated by comparing predicted values of punch force with ex- 
perimental results of reference [5]. Next, a parametric study of punch force and thickness 
strain is carried out to show the effects of process variables namely sheet thickness and 
anisotropy. 


1.3 Structure of Present Work 


The thesis is organized as follows. 

In the second chapter, the mathematical modeling of the deep drawing process is 
presented. Th(' l) 0 \mdary condition applied to the given problem is also presented here. 

In the third chapter, the finite element formulation and the numerical scheme are 
presented. 

In the fourth chapter, the validation and the results for a typical case are presented. 
It also includes a discussion on the parametric study. 

Conclusions and suggestions for further work are presented in chapter five. 
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Chapter 2 


Mathematical Modeling of Large 
Deformation Elasto-Plastic Problem 


In this chapter, mathematical model of a static large deformation elasto-plastic 
problem is developed. The description of motion, stress measures, strain measures and 
constitutive laws used in the formulation of the governing equations for static large 
deformation elasto-plastic problems are presented. 


2.1 Updated Lagrangian Formulation 

In the study of the deformation of a body subjected to external loading, often the original 
undoformed a.nd unstressed state of the body is used for the formulation of its equation 
of mot.ion. This is known as Lagrangian formulation. This formulation is convenient for 
small deformation, which is the case in many engineering problems. In such cases, the de- 
formed c;onliguration does not deviate much from the original one and hence the deforma- 
tion c;an bo described by an infinitesimal strain tensor, for which the strain-displacement 
rclatioiLS arc linear. On the other hand, for large deformation problems, one has to use 
a finite strain mofisure, which is expressed by a nonlinear strain-displacement relation. 
Furthoi-mor<i, the equations of motion when expressed in the reference configuration 
df'.poud on file deformation. Hence, for large deformation problems, the Lagrangian for- 
mulation prove.s to be cumbersome with the governing equations being difficult to solve. 
In such ("us(!,s, one .solves the problem using an incremental method known as updated 
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Figure 2.1: Deformation of a domain 


Lagraiigian formulation. In this formulation, it is assumed that the states of stress and 
deformation of the body are known till the current configuration, say at time t. The 
main ohj(X',tive is then to determine the incremental deformation and stresses during the 
time step A/,, i.e from time t to t+At. Here, the current configuration is used as the 
reference c;onfiguration for obtaining the incremental values. Unlike in the Lagrangian 
formulation, a,n incrciinental strain tensor is used. This methodology is particularly use- 
ful for elfisto-plastic materials because the stress-strain relationship in such materials is 
usually expressed in an incremental fashion. 


2.2 Kinematics of Finite Deformation 


Consider a, domain HI at time t with respect to a given global coordinate system Xi as 
shown in tlic figure 2.1. Let this body deform to domain Let the postion vector 

coi'responding to any particle P be *Xi at time t and at time t + At with relative 

displacem<;nt vector as tAui. The relative displacement gradient tensor at time t -f- At 
is ilu'rcforc defined by 


t+At. p, , 
t ^ V 


d^Xj 


_t+At . 
— t ‘‘'iJ 


( 2 . 1 ) 


11 



Or 


f,+At rp 

t ^ij 


d {hi +t Auj) 
d^Xj 


5ii + 


djAuj 

d^Xi 


( 2 . 2 ) 


When I 'F is non-singular, the polar decomposition theorem [44] allows a decom- 
position of the form: 


'TFij (2.3) 

where is an orthogonal (rotation) tensor representing the material rotation 

and '+^'[7 is a positive definite symmetric tensor called right stretch tensor. The right 
sl,r('t.(il> tensor can be diagonalized by the following transformation to obtain the principal 


stn'tches 



(2.4) 

where ''*'‘^'>,[(^] is orthogonal (transformation) tensor and 


U[i(nosum) 

(2.5) 


The tensor *t^‘R and the transformation matrix are used in the development of 

an objective stress measure while the incremental logarithmic strain tensor is obtained 
from the 


2.3 Stress Measures 


It is esscnitdal that an objective stress measure be used in the incremental theory of 
const, il.utive modeling to account for the rigid body rotation that may accompany defor- 
mation. The Cauchy stress tensor, which has great physical significance, is not objective 
and hence cannot be used directly in a constitutive equation. There are numerous objec- 
tives s(,ress nuiasures (See [45], [46]), each with particular advantages and disadvantages. 
S(jnu! of tlicm arc discussed below: 

One of t,he ccnnmonly used objective stress tensor is the second Piola-Kirchoff stress 
tensor > whicli can be related to the Cauchy stress tensor using the concept 
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of equivalent work between two configurations [47]; 


t.+At a _ P 

‘ “ t^p 




t+At 


^mn 


( 2 . 6 ) 


Here, ‘p and '+^'p are the densities at time t and t+At respectively and t+At'xi,^ denotes 
the derivative of with respect to The second Piola-Kirchoff stress tensor is 

energy conjugate with the Green-Lagrange strain tensor. This energy conjugate pair 
is used in the equation for predicting displacements in a predictor corrector solution 
procedure (see section 3.6.2). 


Allot, hc-r commonly used objective stress rate measure is the Jaumann stress rate 
(T. It is relat.ed to the Cauchy rate a by 


' 4 At = ‘ 4 At aikCWjkAt) a^kCWikAt) (2.7) 

where 

HVijAt = -{tAuij -t Auj^i) ( 2 . 8 ) 

represents tlu' components of the incremental rotation tensor. Here, t+AtAu denotes the 
derivative of tlie incremental displacement vector i^Auwith respect to fX. 

As pointed out by Metzger and Dubey [48], it is important that the stress measure 
used has to he cxnupatible with the constitutive equation in addition to being objective. 
An incremental objective stress measure (as against objective stress rate measures) to be 
usi'd in the generalized Hooke’s law as well as in the elasto-plastic constitutive equation 
is presented Ix’low. 

1'wo Ca,rt,e,sian reference frames (see Fig. 2.2) are used for this purpose and they 

are: 


1. The fixed frame. 

2. A material frame which rotates and translates along with the material particle. 

I’he ma.t,('ria.l frame is defined so as to coincide with the principal axes of the right 
stre(,ch tensor a,t time t SO that the initially orthogonal axes do not get skewed at time 

t -t At . 
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The fixed frame is denoted by subscript F. 

The material frame is denoted by subscript M. 

Hgure 2.2: Fixed and Material Frames of references. 


ddit' first stc'p is the transformation of the components of the Cauchy' stress tensor 
at (.iiiie t from iha fixed frame to the material frame: 

(2.9) 

The incrcMU'iit, in tlie components of the Cauchy stress tensor with respect to the material 
axes, , is added to to obtain the stress components at time t+At with respect 
to the material frame of reference: 

( 2 . 10 ) 

The final trn,nsformation from the frame at time t + At to the fixed frame gives the 
comixments of tTe Caucliy stress tensor at time t + At. 


Tlic ecinatlous (2.9), (2.10) and (2.11) can be combined together to give a relation be- 
tween the et)m])onents of the Cauchy stress tensor at times t and t+At and the increment 
in t.he components with respect to the material frame: 












PJ’ 


( 2 . 12 ) 
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The proposed stress and the logarithmic strain measures satisfy the requirement of objec- 
tivity and lead to a physically consistent application of the usual constitutive equation. 


2.4 Strain Measures 


The incremental small strain is given by 


tAeij = -{tAuij -t- Auj^i) 


(2.13) 


d’lu' increnunit,!!,! Green Lagrange strain is a non-linear function of the displacement and 
is giv(ui by: 


Tlu^se t,wo st.rain measures occur in the virtual work expression at time t + At and its 
transfornm.ti()ii to time t respectively. The components of the Green Lagrange strain 
tcuiKoi- a,r(' inva.riant under rigid body rotation of the material, unlike the small strain 
conipoiunit.s. The following are the disadvantages of using one of the above measures in 
a constitutive law. 

1. The solution obtained is dependent upon the size of the increment in the updated 
Lagrangian formulation unless the increment size is sufficiently small. 

2. The components of the strain tensors do not tend to infinite values when the 
priiuipa,! stretches tend to zero. Therefore, a constitutive law, which ensures that 
t.lu' approjmiate Cauchy stress components tend to negative infinity (as is physically 
rca,listic), even though the strain components remain finite, should be used. This 
difiic.ull.y can be avoided by using a strain measure whose components become 
minus infinity when the principal stretches become zero. 

'rh(i l()garit.hmic strain measure introduced by Dienes [45] is free from the above disad- 
vfuitagcs. d'hc principal incremental logarithmic strain components, which will be used 
in this work, arc; dc.fined by 

tAefj = (2-15) 
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where the are the principal stretches defined by equation (2.5). The logarithmic 

strain has the following additional advantage in elasto-plastic analysis. A loading test 
involving (da.sto-idastic deformation followed by elastic unloading reveals that the slope 
of the ela.sti(, nnlodcling line is the same as that of the initial elastic line only when the 
tiuc sticss and the logarithmic strain measures are used in a constitutive law [49]. 


2.5 Elastic Constitutive Equation 


l he ('la.stic consl.il.ul.ive equation used is the generalized Hooke’s law relating the incre- 
ment in stress components with respect to the material frame and the elastic part 
,As*ij‘ of (.he principal incremental logarithmic strain components: 

= Cfjki tAef. (2.16) 

Tlie l.ensor C''' for the isotropic case is given by 

^ijki ~ X6xj 5]ii -{• l^iSikSji (2.17) 

where A a.nd // iU'c; Ivame’s constants (/r is also called as shear modulus). For anisotropic 
Ixduwior, tlu* tensor C‘‘'‘ has to be evaluated with reference to the directions of the 
principal stretclu's, The stress and strain measures used in a constitutive equation need 
not nec(:sHa,rily l)c energy conjugate with each other. However if it is so, the predicted 
res])ons(^ in a, prc<lic;tor-corrector scheme will be closer to the actual response. 


2.6 ElastoPlastic Constitutive Equation 


As st.rc!ss('s d(’V'(dopcd in a material exceed the yield stress, the elastic constitutive lela- 
tion.ship Ixitwc'cm the st.ress and strain tensors no longer remains valid. A relationship 
bcl.wci'u the stress and strain based on the von Mises yield criterion and isotropic hard- 
(‘niiij’,is pr('scnl.e<l in this section’ . 

' 111 (.lii.s .seel ion the suli.scfiiit/supcr.script denoting time have been omitted for sake of convenience. 
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For an isotropically hardening material, the plastic potential is given by [50] 

= ae,(cry) - ayip). (2.18) 


Note that 

^ = 0 (2.19) 

represents the yield criterion. The plastic potential depends on the Cauchy stress 
tensor aij through it’s second invariant a^g called as equivalent stress and defined by 

(Tij J (2.20) 

wluu'c^ a is the deviatoric part of (Ty. Further, depends on the variable yield stress 
of: the inatc'i'ifd, ay, through a hardening variable p. For the case of strain hardening, p 
is idc!nt,ifi(vl a.s the equivalent plastic strain e^g, and is defined as: 


p s < = / < 


( 2 . 21 ) 


where. 



( 2 . 22 ) 


Here, clef/' is the plastic part of the incremental logarithmic strain tensor defj and the 
integration in eciuation (2,21) is to be carried along the particle path. The dependence 
of ay on p (or el’.g) is normally approximated by a power-law type of relationship: 


= + P-23) 

1 1 ere, {ay),, is the yield stress at zero plastic strain , K is called the hardening coefficient 
and v. is ea,lle(l as I, he hardening exponent. 

The i)liusti(: part of incremental logarithmic strain tensor (de^/) is obtained from 
the plastic pot.ential using the following relation: 

= d\^ ( 2 . 24 ) 

daij 

wlierc dX is a, scahir. This equation is called as the flow rule. Differentiation of equation 
(2.18) with Inspect, to aij gives 
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tlK' definition of (equation 2.22), one can determine dX as 


dX = del^ 


(2.26) 


Fnrtlier, tl.e hardening relationship and the yield condition can be used to express dX 

a,s: 


whorc' 


(l\ 

^ H' ^ H' 


(2.27) 


= ^■'•^4;)”"' (2.28) 

Ls t.lic' slope of t lu' hard('ning curve. Substitution of equations (2.25) and (2.27) in equa- 
tion (2. 2-1) l(>n,ds 1,0 the following constitutive equation: 



SdCTfig ! 


(2.29) 


This constilulive relationship between the deviatoric stress tensor and the plastic part of 
iucrcunental logarithmic strain tensor is not really convenient for the updated Lagrangian 
formulation for which the incremental stress-strain relationship is needed. This can be 
obtained from cxiuation (2.29) as follows: 




da, 


eq 


2H'aeq dan 


da, 


kl- 


Note that, from ccpiations (2.18) and (2.25), one can get 


(2.30) 


^ = -^ = —a'.^. (2.31) 

daki daki 2a eq 

Substitution of equation (2.31) in equation (2.30) leads to the following incremental 



])la.siic stiT'SH sti'aJii relationship; 






kl 


4H'al^ 


da, 


kl- 


(2.32) 


The incienu'iit.ai elastic stress strain relationship (equations (2.16) - (2.17)) can now be 
written as; 


eL 


u,Cy 


-[ 

s'- 


-udakkdij + (1 + i')daij] 


(2.33) 


where de'if is tJu- ('laatic part of deA, E is the Young’s modulus and 
ratio. Adding (lie t.wo relationship equations, one can get 


ly is the Poisson's 



fie)',.' 4 




E 


<%j ^kl 4 - 


1 + t/ 

E 


^ik^jl 4 - 




kl 




da ic 


(2.34) 


This is the iiieix'int'ntal elasto-plastic stress strain relationship needed in the updated 
Lagra,ugia.u ro)'imila.tion. However, it is the following inverse relationship that is more 
useful; 


where’ 




(2.35) 


/ iEP 
^UJkl 


2fi Sik Sji + 


^kt 

1 - 2i/ 


2{3ii + H')alJ- 


(2.36) 


Not.e that, tlu! stress increment appearing in equation (2.35) must be an objective 
st.ix'ss ineixaneiit in the sense that daij must reduce to a zero tensor in the event of the 
iiierciuenl. Ix’iiig a, pure rigid rotation. The incremental objective stress measure to be 
UH('d in the ijreseiit work has been described in section 2.3. 

T'he I'chttiouship (2.35) has been derived assuming the increment size to be small. 
When a, hirgc .size increment is to be used,the relationship (2.35) takes the following form 
wit.li i'<'sp(>('t to the material frame. 
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where 


Aerff = 


r£+At 




(2.37) 




^ji 


+ 


w 


1 -2i/ 


^ij ^kl 


9 mv:,v; 




2{3ij. + Kn{^elg)--^)^a% 


(2.38) 


Here, if ha.s l)een replaced by the expression (2.28) and the left superscript t has been 
added to iiuike it explicit that these quantities are to be evaluated at that time. 


2.7 Incremental Updated Lagrangian Formulation 


The ()hje<;t.iv(' of flic; updated Lagrangian fornaulation is to establish static equilibrium 
in the configuration at time t+At when all static variables at time t are known. The 
principle of virtual work requires that 


where 





R 


(2.39) 


j = Cauchy stress tensor at time t+At, 

M AiV' = volume at time t+At. 

The quantit.y ff' is the virtual linear strain tensor and is defined as 


rt+A(. 


/ d5tAui 


d6tAuj\ 


where 


(2.40) 


” virtual incremental displacement vector at time f, position vector 

a,t tiuKi /. I At. 

Ihirtlua-, ' ' /f is the virtual work of the external forces and is given by 
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<+^<R= [ S (2.41) 

where the traction vector specified on the boundary is given by 

(2.42) 

Here, is the unit outward normal vector to 

The main difficulty in the application of equation (2.39) is that the configuration 
at time t+At i.s unknown. The virtual work expression at time i + At is transformed to 
an integral over the volume at time t [51]. It is assumed that the external load term 
(2.41) is deformation independent for the formulation of the governing equation. The 
expression (2.39) after the transformation becomes 

f '■\‘^'Sij5Ct^'eij)d‘V = R (2.43) 

t\/ 

Where, the virtual Green-Lagrange strain tensor <5 is defined by 

S + t ^Uj,i + t ^Aufej] (2.44) 

The second Piola-Kirchoff stress tensor can be decomposed as 

*TSij =[ Sij +t ASi^ aij +t ASij. (2.45) 

Further, the virtual Green-Lagrange strain tensor can be decomposed as 

= 6{tAeij) = 6{tAeij +t A%-) (2.46) 

where , 

tAsij = -{tAui^j +t Auj^i), (2-47) 

tArjij = -{t^Uk,it^Uk,j)- (2.48) 

Therefore, eciuation (2.43) can be written with incremental decomposition as 
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^SijS{tAeij)d^V + ASij 6{t^r]ij)d^V+ 

{t^ViMV + ^aijd{tAeij)d^V R. (2.49) 

The above equation is simplified by neglecting the second integral, which is a higher 
order term, and approximating tASij as 

'C,%tAekiS{tA6ij)d^V + %j5itArkj)d^V 

+ 1^ 'orij6{tAeij)d^V = (2.50) 

The linearized equation, when solved, will yield only approximate displacement, strain 
and stress fields. The approximate quantities are denoted by a right superscript (1). 
The error duo to the approximation involved is calculated from equation (2.39) as 

Error R - (2-51) 

This error is generally minimized by an iterative predictor-corrector scheme described in 
section 3.6.2. 


2.8 Boundary Conditions for 3-D Deep Drawing of 
Rectangular Sheet 


The sheet or the domain is as shown in figure 2.3. 

The global coordinate system is also shown. Since the problem has two axes of 
symmetry, only one fourth of the sheet (domain) is considered for analysis as shown in 
figure 2.4. 
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Figure 2.3; The complete domain of the Problem. 

2.8.1 Initial Boundary Conditions 

The boundary conditions (both essential and natural) during the starting of the analysis 
are discussed below. 

2. 8. 1.1 The Sheet-Punch Interface 

The sheet-punch interface is represented by face OABC in figure 2.4. The sticking 
friction condition is assumed at the interface. However, a node at the interface may 
loose contact depending upon direction of the vertical reaction force. For the nodes in 
contact, 2 ;-component of incremental displacement vector is specified at the interface. 
For the nodes not in contact, the incremental traction vector is zero. Thus, 

1. < 0(contact nodes). Essential boundary condition: 

tAux = 0 , tAuy = 0, tAuz = specified. (2.52) 

2. > 0 (non contact nodes). Natural boundary condition: 

tAtx = 0 ,tAty = 0, tAt, = 0. (2.53) 
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Figure 2.4: The domain considered for analysis. 


2. 8. 1.2 The Free Surfaces 

The free surfaces, i.e. region where applied traction is zero, are represented by planes 
ABCJIHA,RKLSR,TSLMT,NEDGQPN in figure 2.4. Therefore, the boundary condi- 
tions on these planes are: 

Natural boundary condition: 


tAt^ = 0 = tAt, = 0. 


(2.54) 


2. 8. 1.3 X-Z Plane of Symmetry 

Because of syiuinetry, normal component of incremental displacement vector and shear 
component of incremental traction vector are zero on this boundary. This surface is 
represented by plane RNEDOAHKR in figure 2.4. Thus, 

Natural boundary condition: 

tAtx — 0 , tAtx — 0 , 
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Essential bowtulary condition: 


t^ny — 0 . 


(2.55) 


2. 8. 1.4 Y-Z Plane of Symmetry 

Because of syininetry, normal component of incremental displacement vector and shear 
component ol incremental traction vector are zero on this boundary. This surface is 
represented by plane TMJCODGQT in figure 2.4. Thus, Natural boundary condition: 

,Aty = 0 , tAt, = 0, 

Essential boundary condition: 

tAu^ = 0. (2.56) 


2. 8. 1.5 Sheet-Die Interface 

Let / be the coefficient of friction between the sheet-die interface. Then the boundary 
conditions in .x and y directions are; 

1. No slip: 


iAu^ = 0, tAuy = 0 : if ^Hl +^1 < f \%\ (2.57) 

2. Slip: 

I At, = -fiAt^cose, t^ty = -/tAt^sine : if ^Hl +‘ (2.58) 

where 0 is the angle made by the slip derction with positive x-axis. The boundary 
condition in z-direction is given by: Essential boundary condition: 

tAuz = 0. (2.59) 
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2.8. 1.6 Surface with Blank-Holding Force 

The Blank Holder force is applied on the top portion of the plate represented by the 
planes HIJMLKH in figure 2.4. It is assumed that this force is uniformly distributed 
and it is applied incrementally. Therefore, this boundary condition is 

Natural boundar'y condition: 

tAtx = 0, tAty = 0, tAtz = specified (2.60) 

2.8.2 Incremental Boundary Conditions 

As the deformation progresses in the sheet(domain), the position of the nodes and related 
boundary conditions get affected and therefore an incremental updation of boundary 
condition of nodes in accordance with its position and behavior needs to be done. The 
following section explains how this is achieved. 


2. 8. 2.1 Recognition of State of Nodes 

The whole domain is first divided into regions having similar boundary conditions such 
as those shown in figure 2.5. The nodes are then given different status according to their 
position and region. 

The nod(!s if, at any particular increment, come under the punch radius or die 
radius region arc checked for penetration and are given specified displacement to project 
them onto the [mneh/die profile 


2. 8. 2. 2 Punch Penetration and Specified Displacement 

Suppose a nod<; P, at some increment, penetrates to the position Pi or P 2 (in the corner 
region). The actual location of the point Pi and P 2 are shown in the sectional views 
A-A and B-B. The punch profiles and the position vector (t) of the points Pi and P 2 
with respect to the center of curvature are shown in these sectional views. To avoid the 
penetration, the actual location of the node should be at Pj (or P 2 ) which is obtained 
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Figure 2.5: Subdivision of domain. 


l)y oxlondiiig !_ to the profile. Thus, the boundary condition is applied in the following 
manner: 

• First, the penetration is checked by comparing the magnitude of position vector 
\L\ with t he profile radius. Not e' that t he profile radius is bigger than the punch 
radius in the corner region. 

• The rec|uired specified disi)laeeui(’ut- is then calculated. It- is the difference of the 
position vector of point P (position of the node in the previous increment, not 
shown in figure) and P\ (or P^) 

• The whole incrcuneiit- is t-hc' rc'pi'ated so as t-o gt't- t:he node on t-hc' punch profile. 


2.8.2.3 Die Penetration and Specified Displacement 
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Figure 2.6: Penetration of node in the punch region. 


The penetration of a node in the die region is also checked using a criterion similar to 
l lifU. of .sc'ci iou 2.8.2. 2. Sui)pos(’ a node P peuel rale.s into t.lie die profile to a point Pi or 
P 2 (in the corner region) at some increment. The corresponding projection on the die 
profile along tlic' vector / is shown In- positions P* or P*. 

All t.he nodes lying in t he rc'gion of die raditis are first checked for penetration and 
l.h('n t he requin'd specified disitlaceuK'nl. is calciilat.ed. Th<' re(|uired specified displace- 
ment. is e(|ual t.o t he difference of position vectors of points P and t hat, of point. P* or 
P*. The whoh' increnu'ul. is iIk'h rc'i)eal.('d .so <is t.o gx'l tlu' nock' on the di(' profile. 
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Chapter 3 


Finite Element Formulation 


In this chapter the finite element formulation of the problem is developed. First the 
matrix notation used in the development is discussed in section 3.1. Then the finite 
element eq\ia.tious are derived in section 3.2. ■ The static condensation and other numerical 
aspects of tch solution procedure are explained in the remaining sections. 


3.1 Matrix Notation 

For the 3-D case, the components of the tensors, tAe, tA77 and * Ae-^ are represented 
' in array form a..s follows: 

/{A£:}^ = {t^Exxit ^^yy)t (3-1) 

/{A?]} = Au,x it it ) ‘****}’5 (3.2) 

,{ Ae^}^ = tAe^. jAef;, Ae^y,t As^J. (3.3) 

The components of all stress tensors are written as arrays with respective components 
placed in tin; following manner; 

{o"}^ = \t7xx> ^yyi ^zzi ^yz! ^xz}- (^■^) 

Two matrix forms of the tensor 'Cfjki ( given by equation 2.38) are needed. 
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1. for the constitutive relationship : {{AS"} = j{Ae}, 

2. for the constitutive relationship : t{Ao-^} = ^[C^^]t{A£}^ . 

For isotropic case they are given by: 



In the above ccjuations, 

'^XX ^<^'yy (S'^) 

and a'ij represents the components of the deviatoric part of the Cauchy stress at time t. 
For an isotropic material, [C^] and[C^]' for the 3-D case are given by 

1 — 1 / 2 / 2 / 0 0 0 

2/1 — 2 / 2 / 0 0 0 

2 / 2 / 1 - 2/0 0 0 

0 0 0 0 0 
0 0 0 0 0 
0 0 0 0 0 




1 — 2 / 2 / 2 / 0 0 0 

2/1 — 2 / 2 / 0 0 0 

£ ' u 2/1— 2/0 0 0 

[^'1 ~ (1 4 - 2/)(l - 22/) 0 0 0 1 - 22/' 0 0 

0 0 0 01- 22/0 
0 0 0 0 0 1 - 22 / 


(3.9) 
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Equation (2.50) can be written in the following form owing to the symmetries in 
fAe, {At?, and ‘ cr: 

f (^(aA£}^)‘[C^^]t{Ae}dV+ / 5{t{^rify[T]t{^v}d^y+ 

yf-v Jtv 


f <5(t{A£}^)‘{^T}d‘I/ R. 

Jty 

Here, the ma-trix ' [T] is given by 


where 


1r] = 


‘[E] 0 0 

0 ‘[E] 0 

0 0 ‘[E] 




^ XX 


^xz 

'^xy 

*'^yy 

^(^yz 

^xz 


^CTzz 


(3.10) 


(3.11) 


(3.12) 


3.2 Finite Element Equations 


The domain is discretized into a number of elements and the incremental displacement 
field is approximated over each element by 


AUx 

t{Au} = ^ AUy 

t Arzz 


[$]t{Ati}® 


(3.13) 


where the cknncntal incremental displacement vector t{Au}= for an n-noded element is 
given by 


I { AuY^ ^ {/iAui w AuJ 1 A^l > ’t ’t A«j, ,t Auj | , 


(3.14) 
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The quantities /.AuJ., t^Uy, t^u\ stand for the unknown incremental displacements of 
node i in X',y and z directions respectively and the matrix is defined by 


where 



‘{^ 1 ^ 


(3.15) 


= { tN,, 0, 0, ‘iVa, 0, 0, %, }, 


‘{$2}r={0, ‘ATj, 0, 0, 0, 0,' }, 

W’={0, 0, Wi, 0, 0, 0, }. (3.16) 

The '^Ni, which are functions of (‘x/y/z) are called as shape functions. The number 
of nodes and the corresponding shape functions are chosen on the basis of convergence 
criteria. For the problem under consideration, 8-noded brick element with tri-linear 
shape functions [53] is used. 

The strain field is expressed in terms of the nodal displacements by differentiating 
(3.13) and using the expressions (2.47 and 2.48). This leads to 


,{A£} [5i]i{Aur, 

(3.17) 

t{Av} = [Bjv]t{Au}® 

(3.18) 


where 
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(3A9) 

■{»3};+‘{«i}r, 

'{J-.K+'WS. 

'(■B«r = [ ‘{«>i},.. •{«.}, ., '{* 2 },., ‘{* 2 },, ]■ {3-20) 

Using’ (iciuatious (3.13), (3.17)and (3.18), the contribution to the integral (3.10) 
over a typical clement e ■with volume U® is 

S (.{A«}«^) (jf^^ ‘[S„|’''lrr(BNl/V') ,{Aa}*+ 

i(,{A2.}‘^) (/^_ ■|Bd’'‘[aKt^) =«(.{Au)'') (3.21) 

The contribution to the term is expressed in terms of the elemental external force 
vector using a standard procedure [53]. Since the variation in the displacement 

vector is arbitrary, the above equation can be ■written as 

^[K]%{Aur +Hfr {py ( 3 . 22 ) 

where the elemental stiffness matrix is given by 
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W=nKLr+^[KMLr, (3.23) 

= l^/[BLf'[C^^'f'[BL]d^V^, (3.24) 

'[KnlY=[ '[BNf'[Tf^[BN]d^V^ (3.25) 

and the elemental internal force vector is 

‘{/r=/ W\a]SV\ (3.26) 

Jiye 


The elemental stiffness matrix ‘[ii']® along with the elemental force vectors *{/}® 
and are assembled to obtain the global equation; 


‘[/^]t{Au}+‘{/}=‘+^‘{F}. (3.27) 

Here, is the global stiffness matrix, ‘{/} is the global internal force vector at time 
t and is the global external force vector at time t + At. Decomposing 

equation (3.27) can be written as 

‘[ir],{Au} +‘ {/} =‘ {F} +, {AF}. (3.28) 

Here, ‘{F} is the global external force vector at time t and f{AF} is the global incremental 
force vector from time t to t + At. In updated Lagrangian formulation, it is assumed 
that the equilibrium equations are satisfied exactly at time t. Thus 


‘{/} =‘ {-?’}• (3.29) 

Then, equation (3.28) reduces to 

‘[F]i{Au} =t {AF}. (3.30) 

This equation is solved to obtain the incremental displacement vector t{ Au}. Static 
condensation of. the coefficient matrix and right side vector is carried out before solving 
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this equation. The condensation procedure is explained in the next section. 


3.3 static Condensation Scheme 


Since the contact conditions need to be applied in an iterative manner, solving the 
full storage ecpiations in all the iterations takes an enormously large amount of time, 
especially when the number of degrees of freedom is high. To reduce the solution time, 
it is desirable to condense the coefficient matrix and the right side vector to the size 
involving only those degrees of freedom to which the contact boundary conditions are to 
be applied. Tins condensation procedure is described in the next paragraph. 

First, the total nodes are divided into two categories. The nodes at which contact 
conditions arc to be applied are called the nodes of type 1. The rest of the nodes are 
called as nodcrs of type 2. To facilitate the condensation, the nodes are numbered in 
such a manner that the nodes of type 1 are numbered first and then the nodes of type 2 
are numbered. Then the essential boundary conditions (i.e, the boundary conditions on 
the incremental displacement vector, which is the primary variable in this case) for the 
nodes of type; 2 are applied to equation (3.30). The essential boundary conditions of the 
deep drawing process are described in section 2.9. 

Next, these; equations are partitioned as follows^: 


[Kn] 

[^2l] [K22\ 

where {A'Uj} denotes the incremental displacement vector corresponding to the nodes 
of type 1 and {Aua} contains the incremental displacements of the nodes of type 2. 
Equation (3.31) can be separated as follows: 


{AnJ 1 ^ I {AFi} 
{Au2} j 1 {AF 2 } 


(3.31) 


[i^n]{Aui} + [i^dlAnz} = {AFi}, (3-32) 


[K2i]{Aui} + [ii:22]{Au2} = {AF 2 }. (3-33) 

»In this .section, the subscript/superscript denoting time have been omitted for the sake of convenience 
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Solving for {Av/.^} from equation (3.33), we get 


{Auz} = [K22]-\{AF2} - [i^f2il{Aui}) (3.34) 

Substitution of this expression for {AU 2 } in equation (3.32) and rearrangement of the 
resulting equa.tion leads to the following condensed set of approximate equilibrium equa- 
tions: 


where 


lii1{Aa} = {Af } 


(3.35) 


[K] = [Kn] - [KnM, 

(3.36) 

{AG} = {Aui}, 

(3.37) 

{AF} = {AFi} - [Kn][B], 

(3.38) 

[^] = [iC22]-'[^2l], 

(3.39) 

{5} = [i^22]-'{AT2}. 

(3.40) 


The number of equations in the condensed set (3.35) is equal to the number of 
degrees of frcjedom of type 1. The coefficient matrix [K] and the right side vector {F} 
are evaluated from equations (3.36-3.40) using the partitioned matrices and vectors of 
equation (3.31). Note that while evaluating [^4] and {B}, it is not necessary to invert the 
matrix [^^ 22 ]- Instead, one can solve the following equations by the Gauss elimination 
method: 


[^22] [^] — [-^21]. 


(3.41) 
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(3.42) 


[K22]{B} = {AF 2 }. 

To reduce the computational time further, one can store the upper triangular from of 
[K 22 ]- In this way, one can obtain columns of [.A] by performing the Gauss elimination 
and back substitution operations on the corresponding columns of [K 21 ]. The vector {S} 
is obtained in similar fashion by performing the Gauss elimination and back substitution 
operations on the vector [AF 2 ]. 


3.4 Determination of Stresses 

The evaluation of the stress components (at the Gauss points of the elements) is done 
by the following stepwise procedure. 

1. Calculate the relative deformation gradient using equation (2.2); It should 

be noted that the position vector ‘x corresponds to the equilibrium position. 

2. Decompose the relative deformation gradient as and de- 
termine and using equation (2.4). 

3. Determine the principal incremental logarithmic strain ({Ae-^} using equation 
(2.15). 

4. Calculate t{Ac 7 '^}using equation (2.37). The integration of the constitutive equa- 
tion is performed using the Euler forward technique described below. 

5. Use equation (2.12) to calculate the Cauchy stress components at time t + At. 

3.5 Integration of the Constitutive Equation 

Different teclmiciues exist for the integration of the constitutive equation (2.38). A snnple 
and robust technique is the Euler forward integration scheme described below. 

Assume that the principal incremental strain components have been calculated and 
the state of the Gauss point at time t is known. 
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If the state at time t is elastic, 


1. Calculate the stress increment assuming elastic behavior: 




2. Calculate } using (2.12). 

3. Determine and using (2.20). 

4. If then elastic behavior holds, or if =^<Jeq the Gauss point 

is neutrally loaded. Return. 

5. If >‘ Uy, a transition from elastic to plastic has occurred. Calculate the 


Ratio = 


(Teg) ’ 


change the state to plastic and use the sub-incrementation method. In this method, 
the stress components with respect to the material frame are updated after each 
sub increment by the increment in stress components corresponding to the elasto- 
plastic strain sub increment. The corresponding to the last updated state 

is used: 


(+A( =t+A(. dt {As^} fori = 1, n 

vrheiv. 




£+At ^£ {cT^} + Ratio [C^]^ • 


£+A£ evaluated at {(t^}^°^ 

Use equation (2.11) to find Return. 
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If the state at time t is plastic, the sub incrementation method described in (5) 
above is applied with ratio set to zero. 

This describes the Euler forward integration scheme for the elasto-plastic model 
described in section 2.6. 

3.6 Solution Procedure 

Equation (3.35) represents only an approximate equilibrium equation at time t + At 
(the, approximation is mostly due to the linearization and simplification involved in the 
steps between equations (2.49)and (2.50). A solution of such an approximate equation 
may involve a significant amount of error and, depending on incremental displacement 
step, may become unstable. Therefore, it is necessary to modify equation (3.35) to 
turn it into an iterative problem capable of providing a solution with desirable accuracy. 
Amongst various iterative techniques, the modified Newton-Raphson algorithm is the 
most effective as it offers fast convergence with less computation. This algorithm is 
described in section 3.6.2. 

As stated in section 1.2 of chapter 1, in the present problem, the friction at the 
punch-sheet interface is assumed to be sticking friction. The corresponding contact con- 
dition is applied only in an iterative fashion. Thus, there are two sets of iterations 
in the solution procedure of the present problem. First, for the specified incremental 
force/displacement, the contact iterations are carried out to determine the status (con- 
tact or non-contact) of the interface nodes. This algorithm is described in section 3.6.2. 
Then, corresponding to this status, modified Newton-Raphson iterations are carried out 
to minimize the error arising out of linearization and simplification of the equilibrium 
equation. 

Numerical integration scheme for the elemental coefficient matrices and right side 
vectors and .some divergence handling procedures for the modified Newton-Raphson 
scheme are presented in sections 3.6.3 and 3.6.4 respectively. 
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3.6.1 Contact Algorithm 


This section describes the iterative procedure to apply the contact boundary conditions 
(2.52 and 2.53) to the linearized and condensed set of finite element equations (3.35). 
The iterative algorithm is as follows: 


(I) First iteration: 

1. In the first iteration, it is assumed that all the punch-interface nodes are sticking. 
Thus, the boundary condition given by the equation (2.52) is applied to all nodes. 

2. Next these equations are solved to obtain {tAui}. The vector {tAii 2 } is determined 

from ecination (3.34). Thus, the whole incremental displacement vector is 

known. 

3. Then, the incremental reactions are found by multiplying the original coefficient 
matrix with the incremental displacement vector. 

4. At the end of the contact iteration, the non contact nodes are identified using the 

condition > 0. 

(II) Second iteration: 

1. Now the condition (2.53) is applied to all the non contact nodes. 

2. Then the resulting equations are solved to find {jAui}. The vector {tAu 2 } is found 
as before. 

3. Now the reactions are found by the method described above. 

4. Finally, a check is made to find whether any additional nodes are loosing contact 
using the condition mentioned above. 


(Ill) Further iterations: 

The second step is repeated till there is no change in the contact status between 
two successive iterations. 
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3.6.2 Modified Newton - Raphson Scheme 


The condensed set of approximate equilibrium equations (3.35), when the subscript /superscript 
denoting the time are restored, takes the form: 


^[K]^{Au}=t{AF}. (3.43) 

The modified Newton- Raphson iterative scheme is used to convert this problem into an 
iterative problem to improve the accuracy of the solution. This algorithm can be stated 
as follows. 

1. Form the FEM equation, 

(3.44) 


where 


2. Get the condensed form of the equation 3.44, 

'[i?]dA2}W = 

3. Check for convergence given by the following criterion, 

4. If not converged go to step 1. 

, 5. If converged go to next increiiient. 


(3.45) 


(3.46) 


(3.47) 


The vector {ERR} is called the unbalanced force vector. The right superscript 
on denotes the configuration on which the integration is to be performed for 

the external force vector. For first iterative step, 

>+^<{ERR}^^^ =t {AF} (3.48) 
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3-6.3 Numerical Integration Scheme 


Exact evakuiticn of integrals appearing in element coefficient matrices and right side 
vectors is not always possible because of the complexity of the integrand. In such cases, 
it is natural to seek numerical evaluation of these integral expressions. Numerical inte- 
gration involves approximation of the integrand by a polynomial of appropriate degree, 
because the integral of a polynomial can be evaluated exactly. 

The most commonly used numerical integration method is Gauss-Legendre inte- 
gration method. The expression for the integral over a master element is given by 


Fit V, Od^dvdC = Fit T), Od^dridC 

m n L 

- CkWi>^j‘^k (3.49) 

i=i j=i fc=i 

where rri, n and I denote the number of quadrature points .(Gauss-points) in t V C 
directions re.si>ectively and Ui, ujj and uJk denote the corresponding weights. Two Gauss 
points in each direction are used in this work. 


3.6.4 Divergence Handling Procedures 

The modified Newton-Raphson method diverges in some cases. The following simple but 
fairly effective techniques, for overcoming the divergence, are incorporated in the present 
work. 


1. Line search: When there is slow convergence or divergence in iterative step of an 
increment, the line search is invoked. In this, the right side vector of equation (3.44) 
is scaled by a factor a;. A number of values of on in the range of iamin oimax) 
are used. The value of a/ correponding to the minimum unbalanced force vector 
ERE. is used. Thus, instead of rejecting the whole iterative step, a compromise 
is made to reach the minima for the same iterative step. In the terms of the vector 
'■^^ER.R! for iP'’ increment 

i{ERRY < ai{ERRy : amin <ai< (3-50) 
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2. Incremental Cutting Method: When the Modified Newton Raphson method is 
unable to converge within specified iterative steps, the whole increment is repeated. 
The increment is repeated with the reduced incremental step size. In the terms of 
the specified displacement, this condition can be stated as: 

^^^spedfiedy ^ Oc < 1 (3.51) 

The Inc:rement cutting method is computationally expensive and is used as a last 
resort. 
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Chapter 4 


Results and Discussion 


The finite element model of 3-dimensional deep drawing process developed in the previous 
chapter has been applied to a number of cases involving two materials (Aluminum-killed 
Steel and material in reference [19]) and various sets of process parameters to illustrate 
its apphcability. The domain considered in section (2.8) of chapter 2 is descritized as 
shown in figure (4.1). Anticipating a stress concentration at the corner of the punch and 
in the clearance gap, the descritization is refined in these regions. The validation of the 
analysis is performed in section (4.1). The thickness strain distribution, deformed mesh 
and stress distribution for a typical case is discussed in section (4.2), (4.3) and (4.4). A 
parametric study is done in section (4.5) illustrating the effect of various parameters on 
punch-load and thickness strain. 
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4.1 Validation 


The fem code developed simulating the deep drawing process is validated with respect to 
experimental results given in reference [4]. The material used for this is aluminum-killed 
steel with physical properties and dimensions as follows: 



Figure 4.2: The variation of punch load with punch displacement 


Material: Aluminum killed steel 
Young’s modulus E = 211.744 Gpa, 
Poisson’s ratio v = 0.291, 

Initial Yield Stress {ay)o = 172 Mpa, 
Hardening Coefficient K = 574.59, 
Hardening Exponent n = 0.269; 
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Sheet dimensions (for domain shown in figure (2.4)): 


Blank size s = 55 mm x 55 mm, 

Punch size p = 20mm x 20 mm, 

Die size d = 21.25mm x 21.25 mm, 

Die radius = 5mm, • 

Punch radius Tp = 5mm, 

Corner radius Tc = 3.2mm, 

Sheet thickness t = 0.86mm; 

Blank holding force (for region shown in figure (2.5)) = 500 kg. 


Figure (4.2) shows the variation of the punch load with the punch displacement 
for the above configuration. It can be seen that the punch-load follows closely with the 
experimental result for the initial deformation of the sheet till the punch radius region 
is in fuU contact with the pimch after which there is a deviation in punch force firom 
experimental values with increasing pvmch displacement. The reasons for the deviation 
may be as follows: 

1. As reported in reference [4], at large pimch displacements, the lubricant between 
the sheet and the die may be squeezed out due to which the coefficient of firiction 
may not remain constant as assumed in the analysis. This may further 8dfect the 
stress distribution and the deformation of the domain. 

2. Due to limited computational resource and to limit the computational time, the 
meshing in the blank holder region is not refined. As these elements come under 
the deformation zone, stress developed may be higher than the actual value leading 
to higher punch load. 


.1 


47 



4.2 Thickness Strain Distribution 


Figures (4.3) show the thickness strain distrib\ition in Alnminum killed steel sheet. The 
blank, punch and die sizes as well as blaj}k holding force are the same as in section 4.1. 
In pure deep drawing process, the thickness of the sheet metal cannot change because 
the siirface area does not change [2]. But, in actual practice, the thickness of drawn part 
varies considerably from the initial blank thickness. 



Figure 4.3: Thickness Strain of the sheet metal along the x axis 


The thickness variation depends on the metal, the geometry of the formed part, 
and the drawing technique. The minimum thickness usually occurs within the curved 
portion of the part which connects the bottom with the walls. Figures (4.3) show that, as 
the draw depth is increased, the sheet is most severely stretched in the clearance space. 
The thickness of the sheet under punch is fotind to be less than the original thickness 
whereas, on the die it is slightly greater than the original thickness. 
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4.3 Deformed Mesh 


Figure (4.4) shows the deformed mesh at the end of 14 mm of punch displacement. 



Figure 4.4: The deformed sheet metal 


Aluminum killed steel sheet is used. The blank, punch and die sizes as well as 
blank holding force are the same as in section 4.1. The nodes in the punch radius region 
are considered to be sticking to the punch which is a fair assumption considering the 
large compressive force generated at the sheet-punch interface. 


49 


Figure 4.5: The complete plate showing ear formation at the corner 


The top view for the whole sheet extrapolated from the one fourth portion of the 
sheet is shown in figure (4.5). The formation of ears can be seen at the corner of the 
plate. This phenomenon occurs due to the presence of extra materials. To avoid such 
kind of ear formation, the blank shape needs to be optimized. 
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4.4 Stress Distribution 


Stress Distribution in log scale (N/mm 




Figure 4.6: The Stress distibution in the domain 


The equivalent stress distribution in the sheet during the deformation is shown in 
figure (4.6). It can be seen that the value of equivalent stress is maximum at the punch 
corner and from that point it decreased in all directions to reach a minimum at the sheet 
edges. The stress values are ploted in logarithmic scale. 


4.5 Parametric Study 


Various process parameters affect the final state of deformation and stress in the deep 
drawn cup. These may include sheet thickness, material properties, blank holder force, 
lubrication conditions between the sheet and die etc. The effect of some of these param- 
eters is analyzed in the present work. 
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4.5.1 Effect of Sheet Thickness 

The variation in punch load and thickness strain for different sheet thicknesses is studied 
in this section. The material property used and all other parameters are the same as in 
previous section. The effect of sheet thickness t on punch force is studied by carrying 
out the analysis for five cases (t = 0.50mm, 0.75mm, 0.86mm, 0.95mm and 1.05mm). 
Figure (4.7) shows the variation of punch force with punch displacement for these cases. 
It is observed that the punch force increases with the sheet thickness.Increase in the 
sheet thickness implies increase in the volume. Therefore, a large force is required to 
achieve the same punch displacement as shown in figure (4.7). 



Figure 4.7; The variation of punch load with sheet thickness 


Figure (4.8) illustrates the variation of thickness strain along one of the axis for 
different sheet thicknesses. No uniform variation in thickness is obeserved. However, 
maxim u m thickness strain increases with thickness. 
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Figure 4.8: The variation of thicknes strain with sheet thickness 


4.5.2 Effect of Material Properties 

The punch force and thickness strain distribution are studied for one more material [19] 
the properties of which are as follows: 


Material: 

Young’s modulus E = 71 Gpa, 

Poisson’s ratio u = 0.33, 

Initial Yield Stress (<Ty)o = 132.225 Mpa, 
Hardening Coefficient K — 444.45 Mpa, 
Hardening Exponent n = 0.52836 


The punch force and thickness strain of the sheet for this material are shown in 
figures (4.9) and (4.10). It can be seen that the trends in these figures are the same as 
for earliar material. 
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Figure 4.9: The variation of punch load with material 



Figure 4.10: The variation of thicknes strain with material 
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Chapter 5 


Conclusions and Scope for Future 
Work 


The deep drawing process is one of the most common manufacturing process for sheet 
metal parts. The nature of the deformation involved in this process and the complexity 
in describing the deformation mathematically, renders this process diflScult to analyze 
and simulate. 

A finite element simulation of deep drawing process using solid element is at- 
tempted. The incremental updated Lagrangian formulation is used with new incremental 
objrctive stress measure and logarithmic strain measure to formulate the mathematical 
framework for this simulation. Modified Newton-Raphson method is used to solve the 
resultant non-linear finite element equations. The material is assumed to be elastic- 
plastic strain hardening yielding according to von-Mises criterion. The strain hardening 
behavior is modeled by a power law. The inertial and body forces are neglected. The 
sticking friction is assumed at the punch-sheet interface and sliding Mction is assumed 
at the sheet-die interface. The blank holding force is applied incrementally. 

The validation of the finite element code, thus developed, is done with the available 
experimental results. The parametric study has been carried out to show the effect of 
two process variables namely, sheet thickness and material properties. The following 
conclusions are drawn from the analysis: 

• The punch load it increases with the sheet thickness. 


55 



• Material properties have significant effect on the punch load. Higher punch load is 
required for the material with higher yield stress and larger hardening coefficient. 

• The sheet thickness reduces at the clearance region owing to hi gh tensile stress 
developing in the side wall of the drawn cup. This phenomenon manifests in the 
higher value of the thickness strain, which has a lower value in the clearance region. 

• The equivalent stress developed in the sheet is the largest at the punch corner zone 
and reaches to a minimum value at the free edges of the sheet. 

• For the drawing of square cup drawing from square blank, earring is observed. 

This simulation of deep drawing process is based on many assumptions due to 
which the accurate representation of the process is still incomplete. The present work 
can thus be extended to the following; 

• Most sheet metals exhibit anisotropy in their physical characteristics and it needs 
to be incorporated to describe the stress and strain variation. 

• The present work assumes sticking friction at the punch-sheet interface. SHding 
friction can therefore be included at this zone. Similarly, sUding friction can also 
be incorporated on the curved profiles of the punch and die. 

• The present work can be extended to the deep drawing of arbitrary shape. 

• Failure analysis of deep drawn parts in terms of wrinkling, earing, fracture etc. can 
be studied. 

• Optimization of various process parameters such as blank-holder force, clearance 
gap, blank shape etc. can be included in the scope for further analysis. 
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